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Abstract 

We prove that uniform accuracy of almost second order can be achieved with a finite 
difference method applied to Navier-Stokes flow at low Reynolds number with a moving 
boundary, or interface, creating jumps in the velocity gradient and pressure. Difference 
operators are corrected to 0{h) near the interface using the immersed interface method, 
adding terms related to the jumps, on a regular grid with spacing h and periodic 
boundary conditions. The force at the interface is assumed known within an error 
tolerance; errors in the interface location are not taken into account. The error in 
velocity is shown to be uniformly 0{h‘^\ log/ip), even at grid points near the interface, 
and, up to a constant, the pressure has error 0{h‘^\ log Zip). The proof uses estimates 
for finite difference versions of Poisson and diffusion equations which exhibit a gain in 
regularity in maximum norm. 

1 Introduction 

Recently there has been enormons development in nnmerical methods for flnid flow with 
moving bonndaries or flnid-structnre interaction. Often finite difference methods are nsed on 
a Cartesian grid which does not conform to the moving bonndary. A separate representation 
is used for the boundary, and the effect of the boundary on the fluid must be included. 
For biological models practical applications have most often used the immersed boundary 
method [22] in which the force on the fluid from the boundary is spread to nearby grid points. 
Other methods maintain a sharp interface and are seen to attain about 0{h?) accuracy in 
the velocity. Generally these methods are carefully designed to control the truncation error, 
taking into account the location of the boundary relative to the grid cells. Here we focus 
on the immersed interface method ([131 [IH (13 CS [El [201 [25| |26| [E]) in which difference 
operators for the velocity and pressure are corrected where the stencil crosses the interface 
using jumps in the quantities and their derivatives. Closely related methods use ghost points 
or cut cells. With low to moderate Reynolds number, it is often observed in computations 
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that the velocity error is about 0{h?) even when the truncation error is 0{h) near the 
interface, but this phenomenon has not been explained, and understanding of the solution 
error has come mainly from numerical evidence. 

In this work we estimate the errors in velocity and pressure, uniformly with respect to 
grid points, including those near the moving boundary, in a simple prototype problem using 
the immersed interface method. We neglect possible errors in the boundary location and 
concentrate on the errors in fluid variables brought about by the numerical treatment of the 
force from the moving boundary. Thus for a problem with coupled motion of the fluid and 
moving boundary we are only partially accounting for the errors. We verify analytically the 
principle that 0{h) truncation error at the moving boundary can lead to uniform accuracy 
close to 0{h?). In doing so we elucidate the minimal features necessary to achieve this 
accuracy. This result depends on the effect of diffusion with implicit time stepping and thus 
is signihcant at low Reynolds number. 

We will always measure errors in maximum, or L°°, norm. One reason is that methods in 
use are generally designed to control maximum truncation errors near the moving boundary. 
A second reason is that the errors in the solution are likely to be largest near the boundary, 
and the most meaningful measure of the error is a uniform estimate. We use estimates 
derived in [3l H] for discrete Poisson and diffusion equations with a gain in regularity in 
maximum norm. Although we have chosen to study the immersed interface method, we 
hope that the analytic technique introduced here will be suggestive for the larger class of 
related methods. 

We first state the physical problem. We consider fluid flow in a rectangular region 
[—L, LY in dimension d = 2 or 3, with velocity u and pressure p, both periodic. We suppose 
the moving boundary or interface T is a closed curve in or closed surface in R^. We assume 
the density is constant and the Reynolds number is low to moderate, and for simplicity we 
set both to 1. The fluid flow is determined by the Navier-Stokes equations for the velocity 
and pressure, with a force exerted by the interface on the fluid. 


ut + u- Vu + Vp = Am + /hr , V ■ m = 0 (1.1) 

where / is the force density, hp is the delta function restricting to the surface T, and A = 
is the Laplacian. Equivalently, the equation holds away from T with zero force, and the 
velocity and pressure have the jumps across T 


m] = 0, 

du 

dn 

= 

ftan 1 

1 

III 

e 

d 

(1.2) 

[P] = 

f-n, 

dp 

dn 

^r * ftan 

(1.3) 


Here n is the outward unit normal at T, and [p] = p+ — p- is the difference between the 
outside and inside values at T. (See e.g. [131 US [HI 113 ISSl l26].) 

The operator Vp- is the surface divergence; in R^ it is just the arclength derivative. (E.g, 
see [2], (9.41.1) for the dehnition of Vp- and [26] for a thorough derivation of [dp/dn] .) 
The fact that the pressure p is periodic depends on the facts that 


[u ■ Vm] • n 


0 and 



dS 


0 


(1.4) 
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At each t, the pressure has an indehnite constant; adding po(t) to p{x,t) does not change 
(1.1) or (1.3). 

Typically T moves with the fluid velocity and the force / is determined from the configu¬ 
ration of r, depending on its material properties, e.g. elastic forces, so that / and T depend 
on the fluid variables. In this work we assume the location of T is known exactly and / is 
known within a certain error tolerance; see Theorem 1.1 below. Thus, for the full problem, 
we estimate only the part of the error in velocity and pressure from their direct compu¬ 
tation while neglecting the influence of errors in T. With this qualification, the maximum 
errors in velocity and pressure in the scheme studied here are shown to be 0{h‘^\ loghp) and 
0{h‘^\ logh|^), resp. 

We discretize u, p and their derivatives on a regular grid at points Xj = {ji,j 2 )h or 
js)^ with h = L/N. We use the usual centered differences for the discrete 
gradient V/i, divergence V/^- and Laplacian A/^. All are O(h^) accurate for smooth functions, 
and thus at grid points where the stencil does not cross the interface T. In the immersed 
interface method, the differences are corrected at the irregular points using jumps at T for 
the variables and their partial derivatives. For example, for a function u(x), a; G M, if Xj-i, 
Xj are inside T and Xj+i is outside, T intersects the grid line at x* with Xj < x* < Xj+i, and 

= Xj+i — X*, then 


2h’ ‘ ^ ^ h^[Dv\ + + 0(h^) ( 1 . 5 ) 

+ ” 1+1 _ ^ + '^lD‘v]j + 0(h) ( 1 . 6 ) 

provided v is on each side, with similar but different formulas at Xj+p, e.g., see [TBl 1271 12B] . 
Jumps in the first and second partial derivatives of u and p, needed for corrections here, can 
be found from (1.2), (1.3), as explained in the works cited. 

Next we present the scheme to be analyzed. We choose a time step t = 0{h) and compute 

the velocity u" at time tn = nr. In updating u” to we will assume the force / and 

needed corrections are known up to time tn+i- (If the interface T is updated explicitly, the 
simplest possibility, then F”^^ is found from u”, and it determines the jumps at time tn+i-) 
Assuming is given, we start with 

= —t{u ■ Vm)° — rVp° -I- + tC\ -|- tC-j (1.7) 

Then, with u known up to time u > 1, the new velocity is found from 

We will describe the discretization of each term. Here C\ corrects the approximation 
(^yp+i _ points which are crossed by the interface during the interval tn <t < 

tn+i- This correction is a term proportional to [ui\. Since the velocity is continuous at the 
interface, the material derivative Ut + u-Vu is also continuous, so that [ut\ = -u - [Vu]. (E.g. 
see [13], (33)-(35) or [2S|, Cor. 3.2.) Corrections to other terms at locations crossed by the 
interface during a time interval will be included in C 7 , discussed later. 
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The advection term is extrapolated in time for n > 1, 

(m • ■ VhU^ - ■ VhU^-^ + C 2 (1.9) 

where C 2 is the correction to the centered difference Vh, determined from the jumps [Du] 
and [D'^u], at each time tn and tn-i as in (1.5). Similarly 

with corrections to Ah again determined from [Du], [D'^u] at each time as in (1.6). 

For the pressure term we hrst compute the divergence of (1.9), 

V ■ (m ■ = Vh-{u- + ^4 (1.11) 

where € 4 ^ corrects Vh- for the jumps in [Du], [D'^u]. We then solve the discrete Poisson 
problem 

A^P«+1/2 _ . Vm)’^+V2 _ Q + C's - m (1.12) 

with periodic boundary conditions. Here C 5 corrects Ahp"'~^^^‘^ using the jumps [p], [Dp], [D‘^p] 
and [u ■ Vm]. The last term m is the mean value, or average, of —C 4 + C 5 on the grid. It is 
subtracted so that the right side above has mean value zero and thus is in the range of Ah, 
we will see that the error resulting from m is not signihcant. The periodic solution of (1.12) 
has an indehnite constant; we choose it to have mean value zero. Finally, 

Vp’^+V2 _ + Ce (1.13) 

where Cq corrects using [p], [Dp], [D'^p]. 

For each of the terms (m ■ we compute separately at two time 

levels, including corrections. At a grid point crossed by the interface during the time interval 
we add a correction proportional to the jump in each quantity; see (14),(15) in [13]. These 
terms form Cj. 

We will show that this scheme produces values of u and p that have accuracy slightly less 
than 0{h?) with certain assumptions. We summarize the conclusion as follows. 

Theorem 1.1. Suppose the exact solution of (1.1)-(1.3) is smooth in space-time on each 
side of the interface F for 0 < t < T, and also F is smooth. We neglect any errors in F, 
and we assume that f and Vtanf are known within maximum error Then, with t jh 

constant and h sufficiently small, 

max D) - u^^^^\xj,tn) I < Krh^l log h\^ (1.14) 

j,n 

for tn = nr < T and some constant Kt independent of h. The pressure p” can he found 
as above at time D so that, for some constant pg, depending on h, p^ + pg differs from the 
exact pressure with maximum error hounded by logh|^. 

We have assumed periodic boundary conditions for the computational domain to avoid 
the serious issue of handling the boundary conditions, in order to focus on the accuracy near 
the interface. The difficult question of combining solid boundary conditions with the pressure 
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solution has been dealt with extensively for hnite difference methods including projection 
methods (e.g. [SllTllHllIH]). We expect that in principle the two issues are separate, provided 
the interface is away from the outer boundary. The periodic condition is helpful for the 
analysis in that various operators commute in this case. 

Often a MAC, or staggered, grid is used for velocity and pressure, in the present problem 
( |131 l25] i and others. It allows more natural treatment of computational boundaries. A 
second advantage is that the discrete version of the projection on divergence-free vector 
helds is an exact projection. On the other hand, the simplicity of a single grid, as in the 
present case, is a desirable advantage. We expect that a result similar to the present one 
would hold using a periodic MAC grid. A pressure increment scheme, updating the pressure 
rather than finding it from (1.12), might be used, as in [131 [I3- Such a scheme would be 
different from this one because of the effect of the discrete projection Ei and the present 
analysis would not apply directly. 

For temporal discretization of the diffusion we have chosen to use the Crank-Nicolson 
(CN) method since it is the most familiar second-order accurate method allowing time step 
0 {h) with viscosity and since it is often used with interface simulations and with the pro¬ 
jection method. CN has an important smoothing property in maximum norm, proved in 
[13], related to A-stability. This smoothing is largely responsible for the result stated above. 
However, other methods (called L-stable) such as BDF2 have better smoothing properties. 
The result proved here should also be true for these methods and perhaps can be improved. 

A disadvantage with the use of a single grid is that (V/j-jV^ 7 ^ Ah. Instead (Vh-)'Vh = 
Aq, where Aq is the “wide Laplacian”, a sum of second differences such as 

Vxx{xj) ^ (4h)"^ {vj-2 - 2vj + Vj+ 2 ) (1-15) 

Consequently the discrete version P = I — V/i(A/i)“^V- of the projection onto divergence- 
free vector helds is only an approximate projection, i.e., P'^ ^ P (cf. E |9]). We will see 
that P enters in expressing the error in u. We will hnd that the exact discrete projection 
Po, dehned with (Ao)“^ rather than (A/i)“^, is useful in our estimates. 

In Sec. 2 we collect facts about difference operators in maximum norm relevant to this 
work. We discuss P and Pq. We give an estimate for a discrete Poisson problem with gain 
of regularity. We state estimates for CN time steps, also with gain in regularity. We state a 
lemma which allows a grid function near the interface to be estimated in a lower norm with 
a gain of a factor of h. In Sec. 3 we classify the truncation errors made by the exact solution 
in satisfying the scheme. We can allow errors which, for example, are hrst differences of 
quantities 0{h‘^). In Sec. 4 we estimate the growth of errors in the velocity using the results 
of Secs. 2 and 3 and verify the conclusion above. Finally in the Appendix we give a criterion 
for boundedness of a discrete linear operator in maximum norm and show that a certain 
operator relating P and Pq is bounded. 

Previous analysis of Cartesian grid methods with interfaces has concerned elliptic prob¬ 
lems (Eiisillg]) and linear diffusion (E, Sec. 8 ). The effect of discrete projections on 
accuracy with irregular boundaries was studied in [10]. The case of boundary conditions 
at irregular boundaries in Cartesian grids, rather than interfaces, is quite different and has 
a long history; e.g. see [21], Ch. 6. Analysis of hnite difference methods for the Navier- 
Stokes equations usually assumes smooth, rather than piecewise smooth, solutions and uses 

estimates. 
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We use the letter K for generic constants independent of h. D will be any first spatial 
derivative, and Dh any difference operator, not necessarily centered. However, V/i in a 
gradient or divergence will always mean the centered difference. 

2 Operators on periodic grids 

We collect facts about several difference operators on periodic grid functions that will be 
used in the following arguments. Let Qh be the set of grid points xj = jh in where j is 
a d-tuple of integers with \ju\h < L, 1 < z/ < d, and let be the space of periodic grid 

functions on We always use the maximum norm for such functions: For w G 
||m|| = maxa-.g^h Correspondingly, for a linear operator C on or a subspace, 

||£|| = max||£r(;|| for ||r(;|| = 1. Of course we are primarily interested in how these norms 
depend on h. It will be helpful that the difference operators and inverses we deal with 
commute because of the periodicity. 

The Fourier modes 

Ckixj) = , keZ\ -L/h <K< Ljh (2.1) 

form a basis of They are eigenfunctions for A/j and the “wide Laplacian” Aq, 


AhGk = a{kh)ek, 

4 . 2 

a(kh)= 

D=1 

(2.2) 

AoCk = ao{kh)ek , 

Mkh)= L 

U=1 

(2.3) 


Evidently the null space of A/^ consists of constant functions, the multiples of cq; the null 
space A/q of has dimension 2^ and is spanned by with each = Q oi L/h. (Cf. 
[1].) Each operator is invertible on the subspace of J^{LL}/) spanned by the remaining modes, 
which is also the range of the operator. We will call these subspaces and Xq, respectively. 
Note that Xh = {w G w{xj) = 0}, the subspace with mean value zero. Ah is 

invertible on Xh, and we will write (A/j)“^ on Xh- Similarly we have (Ao)“^ on Xq. If Dh is 
any centered hrst difference and w G then DhW G Xq since Dh is zero on A/q- Thus 

{Ao)~^Dh and {Ah)~^Dh are meaningful for any centered Dh- 

Next we discuss the two discrete versions of the projection on divergence-free vector helds. 
The “exact discrete projection”, again using centered differences V/i, is 

Pqv = V- Xh{Ao)~^Vh ■ V ( 2 . 4 ) 

Since Xh ■ Xh = Aq, Xh ■ Fq = 0, and it follows that (/ — Po)Po = 0, or Pq = Pq and 
(/ — Pq)^ = (F — Pq); that is, Pq and (F — Pq) are exact projections on P{flh)- 
The approximate projection P uses Ah rather than Aq, 

Pv = V - Xh{Ah)~^Xh ■ V ( 2 . 5 ) 

To relate the two, we write 

P = I- XhA-^Xh- = Po + Xh{AQ^ - A-^)Xh- ( 2 . 6 ) 
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and 


P-Po = VH{A,,-Ao)A^^A^^Vh- = {AH-Ao)A^\VhA^^Vh-) = A{I - Po) (2.7) 
where A = {Ah — Ao)Aj{^ so that 

P = Po + A{I-Po), PoP = Po, {I - Po)P = A{I - Po) (2.8) 

The following lemma, proved in the Appendix, tells ns that A is bonnded. 

Lemma 2.1. The operator A = {Ah — Aq)AJ{^ on Xh has ||A|| < K, with K independent 
ofh. 

We have estimates for the inverse Laplacians with gain of two discrete derivatives: 
Lemma 2.2. As operators on Xh 

||A-1 <i7o, \\DhA-^ <Ku \\DIA-^\\<K 2 |logh| (2.9) 

where Dh is any first difference operator and D\ is the product of any two, with constants 
independent of h. The same is true for A( 7 ^ on Xq provided the operators Dh are centered 
differences. 

The statement for Ah is proved in [1], Cor. 3.2 and in an eqnivalent form in [3] The case 
of Aq follows by applying the first case to snbgrids with size 2h. The | logh| factor cannot 
be improved; this can be seen by example. We will nse the fact that V/iAg ^ is bonnded. 
This elliptic estimate applies directly to the projections P and Po, since we can write 


(/ - P)v = Vh{Ah)-^Vh ■ V = Vh{Vh ■ A^i)(n - {v)) (2.10) 

where (v) is the average of v, and similarly for Pq. Then from the lemma we have 

||P|| <77|logh|, llPoll < A^lloghl ( 2 . 11 ) 

We will nse estimates for the resolvent R of Ah and the nth power S'^ of the time-stepping 
operator S for Crank-Nicolson. We dehne 

R={I- fA,)-i , S={I + ^Ah){I - fA,)-' (2.12) 

The following is proved in |1]; see (4.12), (4.13), (7.1), (7.2). 

Lemma 2.3. As operators on P{VLh), 

\\R\\<K,, \\DhR\\ < K2T-P\ (2.13) 

||5"|| < As , \\DhS^R\\ < (2.14) 


where Dh is any first difference and the constants are independent ofh and t. 
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We will need to know that a grid function on the irregular points near the interface is 
almost the discrete divergence of a function which is smaller by a factor of h. The following 
is proved as Lemma 8.1 in [4]. Related statements are Lemma 2.2 in Lemma 2.10 in 
m. and Lemmas 2.2 and 2.6 in [3]. 

Lemma 2.4. Let Th = {tn = nr ,0 < tn < T} and assume r/h is constant. Let he a 
subset of flh X Th such that each {xj,tn) G Xh is within 0{h) ofV{tn)- Let tp he a periodic 
function on Qh x Th which is zero outside ofXh- Then there are periodic grid functions 
0 < z/ < d so that 

d 

(p = $0 + and ||<1),^|| < iLh||(p|| (2.15) 

u=l 

with some constant K independent of h, where D~ is the backward difference in direction v, 
and the norms are the maximum over t^ <T as well as jh G klh- 

Lemma 2.4 applies directly to the correction terms and the remaining truncation errors 
near the interface, since they occur only at grid points withing 0{h) of L. If, for example, 
(p is function of 0 {h) on the irregular points, we express the conclusion briefly as (p = 
0{h?) + DhO{h?) or (p = (J + Dh)0{h?). We might combine this with Lemma 2.3 and use 
the fact that the operators commute to conclude that \\S'^Rp>{-,tk)\\ is 

3 Consistency Estimates 

We estimate the error the exact solution makes in satisfying the scheme (1.7)-(1.13). We 
will denote the exact velocity and pressure as v and q to distinguish them from the computed 
quantities u and p. We will verify that the error has the form 0{h?\ log h\) + DhO{h‘^\ log h|); 
the second term means a difference operator applied to a grid function with the specihed 
bound. 

Lemma 3.1. For the exact velocity v and pressure q we have for n>l 

where the quantities on the right are computed as in (1.9)-(1.13) fromv at times tn-iffnffn+i 
with correction terms. For n = 0, 

v^-v^ = -t{v ■ Vvff - rVq^ + + rCi + rC '7 + (3.2) 

Here = Ef; + DhE^ with = 0{hf\ logh|) for n > 1, fc = 0,1, and = 0{h\ logh|). 

If V and q were smooth across the interface, the scheme would have 0{h?) truncation 
error. To verify the statement, we will consider the corrections at the irregular points near 
the interface and the remaining errors. Since corrections are made to hrst order accuracy, a 
typical truncation error p on the whole grid will have the form 

_ j 0 {h) at an irregular point . ^ 

^ \ 0{h^) at a regular point 


8 


Because of Lemma 2.4, we then have rj = 0{h?) + DhO{h‘^). 

To estimate the errors, we first assume the force / is known exactly, so that the cor¬ 
responding jumps in v, q, and their derivatives are exact. Later we consider the effect 
of errors in /. We will use the superscript ex to denote exact quantities at time tn+ 1 / 2 , 
to distinguish them from quantities computed in the scheme from the exact v. The error 

has the form (3.3) since Ci corrects for [nj if the interface crosses 
the grid point, leaving a remainder at such a point of 0(r) = 0{h). In dealing with other 
terms we will hrst neglect such crossings and return to them afterward. 

The error in the advection term, 

• Vhv^ - ■ Vhv^-^ + C 2 -{v- = 0{h^) (3.4) 

uniformly, since at an irregular point the correction C 2 uses [Dv], [D'^v], leaving remainder 
0{h^/h). Similarly, the error |(A/in" -|- A/in"+^) + — (An)®* has the form (3.3) since the 

remaining error at an irregular grid point after correcting with is 0 {h^/h^). 

Because of (3.4), the error in the divergence 

e^ = Vh-{v- + ^4 - V ■ (n ■ Vn)®* (3.5) 

has the form (3.3), with the correction C 4 similar to C 3 . For the exact pressure g®* we have 

Ag®* = —V ■ (n ■ Vn)®*, with the prescribed jumps in p and dp/dn, so that 

A;,g®* = -V ■ (n ■ Vn)®* + ^5 + £5 (3.6) 

where C 5 corrects the discrete Laplacian with [p], [Dp], [D'^p], [u ■ Vn], and £5 has the form 
(3.3). Combining (3.5) and (3.6) we have 

A,g®* = -V, • (n ■ -C^ + C. + e. + e^ (3.7) 

The pressure corresponding to the scheme is the solution with mean value zero of the 
similar equation 

Anq'^ = -Vh ■ {v • V/,n )”+'/2 _ c-, + ^5 - m (3.8) 

where m is the average of —C 4 + over the periodic box. We check that m is 0{h‘^): In 

(3.7) the averages of A/^g®* and the term are zero, since they are differences. Because 
£4 -h £5 has the form (3.3), and the number of irregular points is 0{h^~'^), its average is 
0{h^), and therefore the same is true for —C 4 -|- C 5 . Now the right sides of (3.7),(3. 8 ) differ 
by an error of the form (3.3), and it follows from Lemma 2.2 that g^ — (g®* — gQ*) = 0{h‘^) 
nniformly, where gg* is the mean value of g®* on the grid. Then V^g^ — V^g*^^ = DhO{h?). 
Also Vhq'^^ + Cq — Vg®* = 0{h) at irregular points, and thus is of the form (3.3). With 
yg»^+i /2 _ v^g^ -f Cg, we combine the last two estimates to conclude that — Vg®* 

is also of the form {Dh + l) 0 {h?). 

For each of the terms [v ■ at a grid point crossed by the 

interface during the time interval, there is an additional correction with the jump in the 
quantity. The remaining error in time discretization at such a point is 0(r) = 0{h), and 
this error again has the form (3.3). 

We now consider the effect of errors in the force /. Suppose at first we assume only that 
the error in / is 0{h?). Then the errors in [Du], ]D'^u], determined from (1.2), are at most 
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0{h?) and 0{h). Then, for example, the error in C 3 is 0{h~^ ■ h ■ h^) = 0{h), which is again 
of the form (3.3) since it occnrs only at irregnlar points. Similarly the error in correcting 
A/jM where the interface crosses is 0{h). The error in Ci or C 2 is 0{h^), and in C 4 it is 0{h). 
These are no larger than corresponding errors already estimated. 

It seems that for the pressnre we need the fnrther assnmption that the error in Vtanf 
is 0{h‘^). The errors in [p], [Dp], [D'^p] are then 0{h‘^), 0{h?), and 0{h), respectively. The 
errors in [Dp], [D'^p] contribnte an error to which is 0{h). The error in [p] contribntes 
an error which conld be 0(1). However, the form of Ahp (1.6) is snch that the correction 
from [p] enters O 5 as a difference, and conseqnently this error in O 5 has the form Dhb for 
b = 0{h). (The difference strnctnre for this correction was noted in [E].) This term b at 
irregnlar points has the form {D^ + l)0{h‘^), according to Lemma 2.4, and thus the error in 
O 5 and Ahp is {D\ + Dh)0{h?). By Lemma 2.2 the resulting error in p, and the contribution 
to — g®* above, is 0(h^| logh|), slightly worse than before. This error term leads to the 
log factor in the statement of the lemma. 

For n = 0 the accuracy of the extrapolation in nV •n +Vg is only 0{h). Since DhO{h^) = 
0 (h), is at most 0 (h\\ogh\). 

4 The error in the solution 

We estimate the growth of the error in the velocity. We set w” = — v^. Rather than 

estimate the maximum of w directly, it seems better to estimate separately g"" = PqW^ and 
= (I — Pq)w'^. Of course ||w”|| < || 2 /"'|| + lk”||, but we cannot bound ||g"'|| by liter'll because 
of the logh factor in the bound (2.11) for the projection. We will see that Pq is useful in 
estimating the nonlinear terms. 

We begin by subtracting the equations (1.8) and (3.1) for and The corrections 
Cl, O 2 , O 3 , Cq, C^ cancel, and the the advection terms give us 

g = ( u - V„n)"+'/^ - (n ■ (4.1) 

where we now use the superscript n + 1/2 as a shorthand for the extrapolation in (1.9). The 
pressures p and g^ are dehned by the similar equations (1.12), (3.8), so that Ah(p — q^) = 
-Vh-g and the gradient term in the equation is 

^ _v,A^V, ■ = -(/ - (4.2) 

We can then combine terms to get 

(u ■ - (v ■ Vhn)’"+^/2 + = Pg ^+^0 ( 4 . 3 ) 

and thus the equation for w has the simple form 

-w^ = + (r/2)(A;,n;”+^ + Ahw^) + re^ (4.4) 

where e” is the error in the n-equation (3.1). We introduce the operators R and S from 
(2.12) and rewrite (4.4) as 

= Sw^ - TRpg^^^!'^ + tRe^ (4.5) 
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For n = 0, = 0 and = 0 since vP = and 


= tRe^ 

(4.6) 

We obtain separate equations for y = Pqw and 2 ; = {I — Po)w by applying Pq and {I — Pq) 
through the tc-equation and using the identities (2.8) for PqP and (/ — Po)P. We get 

yn+l ^ _ ^BP^gn+Xl2 ^ ^p^^E^ 

(4.7) 

and 

^n+l ^ _ p^)gn+V2 ^ _ p^^p^n 

(4.8) 

and for n = 0 

y^ = tPoRe^ , = r(/ - Po)Pe° 

(4.9) 

To estimate the growth of the error, we dehne 


5” = max {\\y^\\ + ||;^-||) , n > 1 

l<m<n 

(4.10) 

We will prove by induction that 


{log h)\ nT<T 

(4.11) 


for some constant Kt and for h sufficiently small. We will then have verihed the estimate 
(1.14) for the error w in velocity. Once the estimate is proved for some n, it follows that 
< 1 for m < n and for sufficiently small h; we will use this below for the nonlinear 
term in the error. 

To estimate we write g = gi + 92 + 93 with 


gi=v-VhW, g2 = wVhV, g3 = w-VhW ( 4 . 12 ) 

It will be important that D^v is uniformly bounded for any hrst difference since v is 
continuous at the interface. We will use the notation B{w) for any bounded linear operator 
applied to w, that is, ||i?(tc)|| < i^||tc|| with constant K independent of h. Thus, for example, 
we can write the difference of a product VjWi as 

Dh{vjWi) = VjDhWi + B{w) ( 4 . 13 ) 

since DhV is bounded. Then for gi and 92 we have 

gn+ 1/2 ^ + B{w^) + DhB{w^-^) + B{w^-^) , = B{w^) + 5 ( m ;”-^) ( 4 . 14 ) 

For the nonlinear term, we can assume by induction, as remarked above, that = 

0(1), and the same for so that H^fsH < iF(||t(;"'|| + ||t(;”“^||). 

The main difficulty in estimating is the effect of Pq on 9 i since Pq has norm 0(| logh|). 
Applying it directly would lose stability. Since we have already estimated g, it is equivalent 
to estimate P^g or (/ — Po)g^ and we choose the latter, (/ — PQ)g = V/iAq ■ g. We note, 
using Lemma 2 . 2 , that V^Aq ^ is a bounded operator, since is a centered difference, and 
for some terms in g it will be enough to estimate the divergence. 
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With g as in (4.12), we begin with gi. Writing (/ — Po)gi = {I — Po){v ■ '^hV + v ■ Vhz), 
we treat the hrst term by using ■ y = 0. We apply Vh- to v ■ VhU obtaining (with sum 
over j, and Dj the centered difference for Xj, 1 < j < d), 

Vh-ivjDjy) = DjVh-ivjy) + Vh-Biy) = Dj{vjVh-y) + iDj + Vh-)Biy) = 0 + DhB{y) (4.15) 
so that by the remark above 


{I-Po){v-VHy) = Df,B{y) 


(4.16) 


For the second term in gi, z is in the range of (/ — Pq) and thus is a discrete gradient, so 
that DjZi = DiZj. Then the divergence is (with sum over i, j) 

Di{vjDjZi) = Di{vjDiZj) = D‘l{vjZj) + DiB{z) = Ao{vjZj) + DiB{z) (4.17) 

Then (J - Pq){v ■ Whz) = W^oivjZj) + DhB{z) = WhivjZj) + DhB{z) = DhB{z). 

The term (/ — Po)g 2 has the form V/iAg ■ {w ■ Vhv) = DhB{w). For g^, again by 
induction VhW = 0(1), and we can treat (/ — Po)g 3 like (/ — Po)g 2 - In summary we have 
shown that 

Po^»^+V2 = (4.18) 

where 

||4>ril <^(lb'"ll + lk”'ll), A; = 0,1; m = n-l,n (4.19) 

and {I — Po)g has the same form. 

We are now ready to prove (4.11) by induction. Since = 0{h\ logh|), it is evident from 
(4.9) that (4.11) holds for n = 1. We assume it is true for n and prove it for n + 1. Here 
and below we use the fact that ||Po|| < K\ logh|. From the y-equation (4.7) we have 


= -r ^ ^ r^F^-^PPge^ = Si + S 2 (4.20) 

£=1 i=0 

For Si we use (4.18),(4.19) and (2.14) to estimate for £ < n — 1 

\S'^-^Dhm[\ < K{{n - ( 4 . 21 ) 

and similarly for other terms, while for i = n we use (2.13) to get 

\Dhm^\ < Kr-^B^n (4 22 ) 


Then 


Sil <7^1 


n—1 

Y,({n - (-)r))-‘/2^V + 

l=\ 



(4.23) 


For S 2 we recall from Sec. 3 that = {I+Dh)0{h?‘ \ log h|) for n > 1, and we again use (2.14) 
for 1 < £ < n — 1 and (2.13) for i = n. We treat i = 0 separately using = 0{h\ logh|) 
and (2.14). Then IS 2 I is bounded by a constant times 


n—1 

^((n — £)r))“^'^^h^(logh)^r + r“^'^^h^(logh)^r + h(logh)^r < Kh^{logh)‘^ (4.24) 
e=i 
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since the sum approximates an integrable function of time. Estimates for are entirely 
similar, since A is bounded independent of h, according to Lemma 2.1, and adding the two 
inequalities gives 


n—1 

^"+1 < + K 2 h^{\oghY (4.25) 

t=i 

To simplify this we use Holder’s inequality to obtain 

( n—1 \ 2 n 

^((n-£)r))-2/V + r-2/Vj '^{5^fT + K'^{h^{\oghff (4.26) 

t=i ' i=\ 

The hrst sum is uniformly bounded, and the inequality has the form 

n 

^n+i < A^k^t + B (4.27) 

i=i 

with = (^"')^ and B = (h^(logh)^)^. It follows easily from this and the estimate for 
that h™ has the bound (4.11) for m < n + 1, provided (n + l)r < T. 

Finally we discuss the accuracy of the pressure. We can compute the pressure p" at time 
tn as in ( 1 . 12 ) using (m • Vhu)"" and the corrections € 4 ^, C 5 . As in Sec. 3, — {q — go) = 

0(h^|logh|), where q is the exact pressure, q^ is computed from the exact velocity v and 
corrections as in (3.8) at time tm and go is the mean value of g on the grid. Combining the 
estimate (4.11) or (1.14) for w with (4.14) for g, we now have g = {Dh + l)0{h‘^\ loghp). 
Lemma 2.2 then gives 

p^-q^ = ^l^Vh-{Dh + I)0{e\\ogh\^) = 0{h^\\ogh\^) (4.28) 

Thus — {q — go) = 0{h‘^\ logh|^), as stated. 

A Appendix 

We prove a simple criterion for the boundedness of Fourier multiplier operators on periodic 
grid functions in maximum norm. We use this statement to prove Lemma 2.1. Related 
statements are given in 0, Ch.l and [ 12 ], Sec. 2.5.1. 

We will assume L = tt for convenience and vr/h is an integer. Let Id be the set of integer 
d-tuples j with \ju\ < 'n'/h, 1 < u < d. Then = hid and J^(fl/i) is the space of periodic 
grid functions. For p e we have the discrete Fourier transform and inverse 

m = , vUh) = (A.l) 

jg/d ^g/d 

and the isometry 

= (27r)-''^|,3(«:)|V. (A.2) 

jg/d fcg/d 
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Lemma A.l. Suppose an operator A is defined on iF{Qh) by 


{Aipy{k) = a{kh)fi{k) 


(A.3) 


where a is a function of f E with period 27i in each 1 < z/ < c?. Let ||A|| be the norm 
of A as an operator on iF{klh) with maximum norm. Then 

/ \ d/is / \ (l-d/2s)/2 

< /<■ ^ {\(D:ra(kh)\^ + \a(khf)h'‘] Y. W(kh)\V) (A.4) 


fce/'Z 




keid- 


where s is an integer, s > d/2, is the forward divided difference operator in direction 

v, a sum over 1 < u < d is implied, and K is independent of h. 

Proof. We can write A as a discrete convolution with the inverse transform Oj of a{kh), 

Aip = aj_iip{ih ), Oj = {2n)~^a{kh)e^^^^h^ (A.5) 

k&i'^ 


so that the operator norm of A on has the bound 


j^jd 


Ajj I . 


(A.6) 


We will temporarily extend this sum to all j G with aj = 0 for j ^ and estimate as in 
the proof of Sobolev’s theorem and particularly as in [5], Ch. 1, Thm. 3.1. With R to be 
chosen we have 


El 


an = 


E Mmn- + E 

\j\>R \j\<R 

s (E E \vm 


\j\>R 

< 


21 4 12s 


1/2 


MEi 

\j\<R 


1/2 


E 


1/2 


x\ dx ] Ml + 


1/2 


lil>R ^ Eil<« ^ ^\j\<R 

r- \ 1/2 

dx] Mo 

^\x\>R—\/d J \J \x\<R-\-y/d / 

< Ki{R- + Ko{R + y/df/^Mo (A.7) 


where 

Ml = E M, Ml = ^ |%Pbf + 1). (A.8) 

We choose R = Vd+ {Mi/Mof/fi so that R > y/d + 1 and {R + y/d)/{R — y/d) is bounded 
above. Then both terms at the end of (A.7) have , and the inequality simplihes 

to 


l®il ^ (1 + b1 

j'eZ'Z 


djAs 

..M / ^ 

j'eZ'Z 


{l-d/2s)/2 


(A.9) 
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Next we relate \ j\^aj to {D^Ya{kh). A summation by parts using periodicity shows that 
(27r)-‘^ hfaj , f3 = - l)/h (A.IO) 

k&I'i 

SO that D'^a{kh) is the transform of the last term, and thus by (A.2) 

|/3(j„ft,fc)|"*|o,f = (211)-“ Y, \(D:)‘a(kh)\'‘h“. (A.11) 

jg/d fcg/d 

We note that ~ 1| = 2| sin(jj,/i/2)| > c\jyh\ since \ivh)\ < it. Thus is 

bounded by the right side of (A. 11 ), and summing over u we get 

d 

< A-5;^IW)V(fcA)|V. (A.12) 

jg/d i/=l fcg/d 

Finally, combining (A. 6 ),(A.9),(A. 12 ) gives the conclusion (A.4). □ 

Proof of Lemma 2.1. For simplicity we assume L = tt. From (2.2), (2.3) the Fourier 
symbol of A in case d = 2 is 

= (4 + sD/isl + .sj) , .Sj = sm{Q/2) (A.13) 

where f = kh, and similarly for d = 3. We set (t(0) = 0 so that A is extended to all of J^(fl/i). 
For |,^j)| < TT, |sj| > c|.^j|, and a and da/dfj are bounded. It is easy to check that d'^a/dfj 
is bounded for 7 ^ 0. It follows that the second difference D'^a is bounded at grid points not 
adjacent to 0. However, for grid points near 0, a = OifY), so that is bounded in that 
case also. Thus the right side of (5.4) is bounded independent of h with s = 2, and Lemma 
A.l ensures that the same is true for ||A||. 
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